Flowpipe construction consists of under or over-approximating the sets of states reachable by dynamical systems. Recently a method has been developed for the class of polynomial ODEs with uncertain initial states (see [1], abbreviated XFZ18under). This method consists of reducing the Hamilton-Jacobi-Bellman equation to a hierarchy of semidefinite programs.
In this notebook we consider the problem of approximating the flowpipe of a system of polynomial ODEs using XFZ18under. This is a Julia implementation that relies on the JuMP ecosystem (JuMP, PolyJuMP, SumOfSquares, MathProgInterface, MathOptInterfaceMosek) and the JuliaAlgebra ecosystem (MultivariatePolynomials, DynamicPolynomials). The implementation is evaluated on a set of standard benchmarks from formal verification and control engineering domains.
References:
Quick links to documentation:
Consider the following van-der-Pol system:
$$ \dot{x}_1 = x_2 \\ \dot{x}_2 = -0.2x_1 + x_2 - 0.2x_1^2 x_2 $$
using MultivariatePolynomials,
JuMP,
PolyJuMP,
SumOfSquares,
DynamicPolynomials,
MathOptInterfaceMosek,
MathematicalSystems
const ∂ = differentiate
# symbolic variables
@polyvar x₁ x₂ t
# time duration (scaled, see dynamics below)
T = 1.0
# dynamics
f = 2 * [x₂, -0.2*x₁ + x₂ - 0.2*x₁^2*x₂]
# set of initial states X₀ = {x: V₀(x) <= 0}
V₀ = x₁^2 + x₂^2 - 0.25
# constraints Y = {x: g(x) >= 0} compact search space Y x [0, T]
g = 9 - x₁^2 - x₂^2
# degree of the relaxation
k = 12
# monomial vector up to order k, 0 <= sum_i alpha_i <= k, if alpha_i is the exponent of x_i
X = monomials([x₁, x₂], 0:k)
XT = monomials([x₁, x₂, t], 0:k)
# create a SOS JuMP model to solve with Mosek
model = SOSModel(with_optimizer(MosekOptimizer))
# add unknown Φ to the model
@variable(model, Φ, Poly(XT))
# jacobian
∂t = α -> ∂(α, t)
∂xf = α -> ∂(α, x₁) * f[1] + ∂(α, x₂) * f[2]
LΦ = ∂t(Φ) + ∂xf(Φ)
# Φ(x, t) at time 0
Φ₀ = subs(Φ, t => 0.)
# scalar variable
@variable(model, ϵ)
dom1 = @set t*(T-t) >= 0 && g >= 0
dom2 = @set g >= 0
@constraint(model, ϵ >= 0.)
@constraint(model, LΦ ∈ SOSCone(), domain = dom1)
@constraint(model, ϵ - LΦ ∈ SOSCone(), domain = dom1)
@constraint(model, Φ₀ - V₀ ∈ SOSCone(), domain = dom2)
@constraint(model, ϵ + V₀ - Φ₀ ∈ SOSCone(), domain = dom2)
@objective(model, Min, ϵ)
optimize!(model)
println("Relaxation order : k = $k")
println("JuMP.termination_status(model) = ", JuMP.termination_status(model))
println("JuMP.primal_status(model) = ", JuMP.primal_status(model))
println("JuMP.dual_status(model) = ", JuMP.dual_status(model))
println("JuMP.objective_bound(model) = ", JuMP.objective_bound(model))
println("JuMP.objective_value(model) = ", JuMP.objective_value(model))
# Recovering the solution:
ϵopt = JuMP.objective_value(model)
# Punder <= 0
Punder = subs(JuMP.value(model[:Φ]), t => T)
# Pover <= 0
Pover = subs(JuMP.value(model[:Φ]), t => T) - ϵopt * (T+1)
using ImplicitEquations, Plots
gr() # or pyplot()
G = plot()
_Punder(x, y) = sum([Punder.a[i]*x^exponents(mi)[1]* y^exponents(mi)[2] for (i, mi) in enumerate(Punder.x)])
plot!(G, _Punder ⩵ 0., xlims=(-10, 10), ylims=(-10, 10), color="red")
_Pover(x, y) = sum([Pover.a[i]*x^exponents(mi)[1] * y^exponents(mi)[2] for (i, mi) in enumerate(Pover.x)])
plot!(G, _Pover ⩵ 0., xlims=(-10, 10), ylims=(-10, 10), color="blue")
G
G = plot()
_Punder(x, y) = sum([Punder.a[i]*x^exponents(mi)[1]* y^exponents(mi)[2] for (i, mi) in enumerate(Punder.x)])
Gu = plot(_Punder ≪ 0., xlims=(-8, 8), ylims=(-8, 8))
_Pover(x, y) = sum([Pover.a[i]*x^exponents(mi)[1] * y^exponents(mi)[2] for (i, mi) in enumerate(Pover.x)])
Go = plot(_Pover ≪ 0, xlims=(-8, 8), ylims=(-8, 8))
plot(Gu, Go)
| Package | k | Constraints | Scalar variables | Matrix variables | Time(s) |
|---|---|---|---|---|---|
| JuMP | 2 | 245 | 175 | 8 | < 1 |
| YALMIP | 2 | 152 | 63 | 10 | < 1 |
| JuMP | 3 | 893 | 715 | 10 | ~ 1 |
| YALMIP | 3 | 254 | 121 | 10 | ~ 1 |
| JuMP | 4 | 893 | 730 | 10 | ~1 |
| YALMIP | 4 | 394 | 206 | 10 | 1.18 |
| JuMP | 5 | 2639 | 2309 | 10 | 1.17 |
| YALMIP | 5 | 578 | 323 | 10 | 0.11 |
| JuMP | 6 | 2639 | 2337 | 10 | 0.90 |
| YALMIP | 6 | 812 | 477 | 10 | 1.10 |
| JuMP | 7 | 6725 | 6183 | 10 | 5.62 |
| YALMIP | 7 | 1102 | 673 | 10 | 1.52 |
| JuMP | 8 | 6725 | 6228 | 10 | 6.06 |
| YALMIP | 8 | 1454 | 916 | 10 | 1.10 |
| JuMP | 9 | 15269 | 14447 | 10 | 41.2 |
| YALMIP | 9 | 1874 | 1211 | 10 | 1.58 |
| JuMP | 10 | 15269 | 14513 | 10 | 38.1 |
| YALMIP | 10 | 2368 | 1563 | 10 | 2.73 |
| JuMP | 11 | 31617 | 30439 | 10 | 244 |
| YALMIP | 11 | 2942 | 1977 | 10 | 2.30 |
| JuMP | 12 | 31617 | 30530 | 10 | 231 |
| YALMIP | 12 | 3602 | 2458 | 10 | 6.57 |